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Abstract 

In spite of the recent interest and advances in linear controllability of complex networks, controlling 
nonlinear network dynamics remains to be an outstanding problem, we develop an experimentally 
feasible control framework for nonlinear dynamical networks that exhibit multistability (multiple co¬ 
existing final states or attractors), which are representative of, e.g., gene regulatory networks (GRNs). 
The control objective is to apply parameter perturbation to drive the system from one attractor to 
another, assuming that the former is undesired and the latter is desired. To make our framework 
practically useful, we consider restricted parameter perturbation by imposing the following two con¬ 
straints: (a) it must be experimentally realizable and (b) it is applied only temporarily. We introduce 
the concept of attractor network, in which the nodes are the distinct attractors of the system, and there 
is a directional link from one attractor to another if the system can be driven from the former to the 
latter using restricted control perturbation. Introduction of the attractor network allows us to formu¬ 
late a controllability framework for nonlinear dynamical networks: a network is more controllable 
if the underlying attractor network is more strongly connected, which can be quantified. We demon¬ 
strate our control framework using examples from various models of experimental GRNs. A finding 
is that, due to nonlinearity, noise can counter-intuitively facilitate control of the network dynamics. 
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An outstanding problem in interdisciplinary research is to control nonlinear dynamics on com¬ 
plex networks. Indeed, the physical world in which we live is nonlinear, and complex networks 
are ubiquitous in a variety of natural, social, economical, and man-made systems. Dynamical 
processes on complex networks are thus expected to be generically nonlinear. While the ultimate 
goal to study complex systems is to control them, the coupling between nonlinear dynamics and 
complex network structures presents tremendous challenges to our ability to formulate effective 
control methodologies. In spite of the rapid development of network science and engineering to¬ 
ward understanding, analyzing and predicting the dynamics of large complex network systems 
in the past fifteen years, the problem of controlling nonlinear dynamical networks has remained 
open. 

In the past several years, a framework for determining the linear controllability of network 
based on traditional control and graph theories emerged iffl, leading to a quantitative un¬ 
derstanding of the effect of network structure on its controllability. In particular, a structural- 
controllability framework was proposed [0], revealing that the ability to steer a complex network 
toward any desired state, as measured by the minimum number of driver nodes, is determined 
by the set of maximum matching, which is the maximum set of links that do not share starting 
or ending nodes. A main result was that the number of driver nodes required for full control is 
determined by the network’s degree distribution [0]. The framework was established for weighted 
and directed networks. An alternative framework, the exact-controllability framework, was subse¬ 
quently formulated |0], which was based on the principle of maximum multiplicity to identify the 
minimum set of driver nodes required to achieve full control of networks with arbitrary structures 
and link-weight distributions. Generally, the deficiency of such rigorous mathematical frameworks 
of controllability is that the nodal dynamical processes must be assumed to be linear. For nonlinear 
nodal dynamics, the mathematical framework on which the controllability theories based, namely 
the classic Kalman’s controllability rank condition liMi, is not applicable. At the present there 
is no known theoretical framework for controlling nonlinear dynamics on complex networks. 

Due to the high dimensionality of nonlinear dynamical networks and the rich variety of behav¬ 
iors that they can exhibit, it may be prohibitively difficult to develop a control framework that is 
universally applicable to different kinds of network dynamics. In particular, the classic definition 
of linear controllability, i.e., a network system is controllable if it can be driven from an arbitrary 
initial state to an arbitrary final state in finite time, is generally not applicable to nonlinear dy¬ 
namical networks. Instead, controlling collective dynamical behaviors may be more pertinent and 
realistic [2^23]. Our viewpoint is that, for nonlinear dynamical networks, control strategies may 
need to be specific and system dependent. The purpose of this paper is to articulate control strate¬ 
gies and develop controllability framework for nonlinear networks that exhibit multistability. A 
defining characteristic of such systems is that, for a realistic parameter setting, there are multiple 
coexisting attractors in the phase space ll24l - l29ll . The goal is to drive the system from one attractor 
to another using physically meaningful, temporary and finite parameter perturbations, assuming 
that the system is likely to evolve into an undesired state (attractor) or the system is already in 
such a state, and one wishes to implement control to bring the system out of the undesired state 
and steer it into a desired one. We note that dynamical systems with multistability are ubiquitous 


in the real world ranging from biological and ecological to physical systems yC 

In biology, nonlinear dynamical networks with multiple attractors have been employed to un¬ 
derstand fundamental phenomena such as cancer mechanisms BTII . cell fate differentiation [138 - 


4lh . and cell cycle control 11421 l43h . For example, boolean network models were used to study the 
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gene evolution [|44|] . attractor number variation with asynchronous stochastic updati^ 114511 . gene 
expression in the state space and organism system grow th rate improvement iSi. Another 
approach is to abstract key regulation genetic networks 11471. 14811 (or motifs) from all associated 
interactions, and to employ synthetic biology to modify, control and finally understand the bi¬ 
ological mechanisms within these complicated systems 113 SL 14211 . An earlier application of this 
approach led to a good understanding of the ubiquitous phenomenon of bistability in biological 
systems ll^ . where there are typically limit cycle attractors and, during cell cycle control, noise 
can trigger a differentiation process by driving the system from a limit circle to another steady 


state attractor 03 811 . Generally speakii^ there are two candidate mechanisms for transition or 
switching between different attractors ll40ll : through signals transmitted between cells and through 
noise, which were demonstrated recently using synthetic genetic circuits [50, U]. More recently, a 


detailed numerical study was carried out of how signal-induced bifurcations in a tri-stable genetic 
circuit can lead to transitions among different cell types [14ill . 

In this paper, we develop a controllability framework for nonlinear dynamical networks based 
on the concept of attractor networks [l52ll . An attractor network is defined in the phase space of the 
underlying nonlinear system, in which each node represents an attractor and a directed edge from 
one node to another indicates that the system can be driven from the former to the latter using 
experimentally feasible, temporary, and finite parameter changes. A well connected attractor 
network implies a strong feasibility that the system can be controlled to reach a desired attractor. 
The connectivity of the attractor network can then be used to characterize the controllability of 
the nonlinear network. More specifically, for a given pair of attractors, the weighted shortest path 
between them in the attractor network is an indicator of the physical feasibility of the associated 
transition. We use gene regulatory networks (GRNs) to demonstrate our control framework, 
which includes low-dimensional, experimentally realizable synthetic gene circuits and a realistic 
T-cell cancer network of 60 nodes. A finding is that noise can counter-intuitively enhance the 
controllability of a nonlinear dynamical network. We emphasize that the development of our 
nonlinear control framework is based entirely on physical considerations, rendering feasible 
experimental verification. 


Results 


A complex, nonlinear dynamical network of N variables can be described by a set of N 
coupled differential equations: 

x = F(x,/i), (1) 

where x G denotes the A^-dimensional state variable, F(x,/r) is the nonlinear vector field, 
and yU G represents the set of coupling parameters. In a GRN, the nodal dynamics is typically 
one dimensional. For simplicity, we assume that this is the case to be treated so that the size 
of the network represented by Eq. ([T]) is N. From consideration of realistic GRNs, we assume 
that the coupling parameters can be adjusted externally, which are effectively the set of control 
parameters. Specifically, in a GRN, the various coupling strengths among the nodes (genes) can 
be experimentally and systematically varied through application of specific targeted drugs. At a 
larger scale, the fate of a cell can be controlled by adding drugs to the cell-growth environment, 
which adjust the interaction parameters in the underlying network ll^ . While dynamical variables 
themselves can also be perturbed for the purpose of control, for GRNs this is unrealistic. For this 
reason the scenario of perturbing dynamical variables will not be considered in this paper. 
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FIG. 1: T-LGL survival signaling system and its attractor network, (a) Structure of T-LGL signaling 
network: each node is labeled with its generic name, and the arrowhead and diamond-head edges represent 
excitation and inhibition regulations, respectively. The inhibitory edges from “Apoptosis” to other nodes 
are not shown (for clarity), (b) Attractor network of the T-cell network, which contains three nodes: two 
cancerous states denoted as Ci and C 2 and a normal state denoted as N. Our detailed computations reveal 
that parameter perturbations on 48 edges can drive the system from a cancerous state to the normal state, 
which are indicated with dark dashed lines, whereas the remaining edges in the network are specified with 
light solid lines. The two directed edges in the attractor network are multiple, containing altogether 48 
individual edges corresponding to controlling the 48 solid-line edges in the original network. 


We focus on nonlinear dynamical networks with multiple coexisting attractors. For a given set 
of parameters p, the multiple attractors (e.g., stable steady states) and the corresponding basins 
are fixed. For a given initial condition, the system will approach one of the attractors. Each 
attractor has specific biological significance, which can be regarded as either desired or undesired, 
depending on the particular function of interest. Suppose, without any control, the system is in an 
undesired attractor or is in its basin of attraction. The question is how to steer the system from 
the undesired state to a desired state by means of temporal and small parameter variations that are 
experimentally feasible. 

Control principle based on bifurcation. To motivate the development of a feasible control 
principle, we consider the simple case where the system is near a bifurcation and control is to be 
applied to drive the system from one attractor to another through temporal perturbation to a single 
parameter. That is, the parameter variation is turned on and takes effect for a finite (typically short) 
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duration of time. After control perturbation is withdrawn, the system is restored to its parameter 
setting before control but its state has been changed: it will be in the desired attractor. Let /tq be 
the initial parameter value and the system is in an undesired attractor denoted as x*, and let 
be the desired attractor that the system is driven to. Imposing control means that we change the 
parameter from hq to yUi. The dynamical mechanism to drive the system out of the initial attractor 
is bifurcations, e.g., a saddle-node bifurcation at which the original attractor disappears and its 
basin is absorbed into that of an intermediate attractor [l43n . denoted as x^. Turning on control 
to change /i from /xq to thus makes the system approach x^. This process continues until the 
system falls into the original basin of x^, at which point the control parameter is reset to its original 
value /io so that the system will approach the desired attractor x^. Success of control relies on the 
existence of a “path” from the initial attractor to the final one through a number of intermediate 
attractors. If a single parameter is unable to establish such a path, variations in multiple parameters 
can be considered, provided that such parameter adjustments are experimentally realizable. For a 
biological network, this can be achieved through application of a combined set of drugs 115 31.15411. 
However, even when potential complications induced by inter-drug interactions are neglected, the 
search space for suitable parameter perturbation can be prohibitively large if we allow all available 
parameters to be adjusted simultaneously. We demonstrate below that this challenge can be met 
by constructing an attractor network for the underlying system. 

Attractor networks: an example of T-cell network For a complex, nonlinear dynamical net¬ 
work, an attractor network can be constructed by defining each of all possible attractors of the 
system as a node. There exists a directed link from one node to another if an experimentally ac¬ 
cessible parameter of the system can be adjusted to drive or control the system from the former 
to the latter. There can be multiple edges from one node to another, if there are multiple parame¬ 
ters, each enabling control. Starting from an initial attractor, one can identify, using all accessible 
parameters with variations in physically reasonable ranges, a set of attractors that the system can 
be driven into. Repeating this procedure for all attractors in the system, we build up an attractor 
network that provides a blueprint for driving the whole networked system from any attractor to any 
other attractor in the system, assuming at the time the latter attractor would lead to desired function 
of the system as a whole. All these can be done using relatively small parameter perturbations. 

To demonstrate the construction of an attractor network, we take as an example a realistic 
biological network, T-cell in large granular lymphocyte leukemia associated with blood cancer. 
Specifically, apoptosis signaling of the T-cell can be described by a network model: T-cell survival 
signaling network ISSioa, which has 60 nodes and 142 regulatory edges, as shown in Fig. [fla). 
Nodes in the network represent proteins and transcripts, and the edges correspond to either activa¬ 
tion or inhibitory regulations. Experimentally, it was found that there are three attractors for this 
biophysically detailed network |55, 5^. Among the three attractors, two correspond to two dis¬ 
tinct cancerous states (denoted as Ci and C 2 ) and one is associated with a normal state (denoted 
as N). By translating the Boolean rules into a continuous form using the method in 11571 15811 and 
setting the strength of each edge to unity, one can obtain a set of nonlinear dynamical equations for 
the entire network system. Direct simulation of the model revealed that there are three stable fixed 
point attractors, in agreement with the experimental observation |55, 5^. The attractor network 
is thus quite simple: it has three nodes only, as shown in Fig.fTJb). Testing all experimentally ad¬ 
justable parameters, we find multiple edges from each cancerous attractor to the normal one (see 
Supplementary Table 1). Since the goal of control is to bring the system from one of the cancerous 
states to the normal one, it is not necessary (or biologically meaningful) to test whether parameter 
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FIG. 2: Relationship between edge control strength and minimal control time. For the T-cell network, 
(a) an inverted rectangular control signal of duration r and amplitude Afj, = — ^o|^ where /io is the 

original parameter value. A saddle-node bifurcation occurs for // = /Tc, so Ae = /Tc — /in is the excess 
amount of the parameter change over the critical value /ic. (b,c) Minimal control time Tm versus /in, where 
parameter control is applied to the activation edge from node “SIP” to node “PDGFR” and to the inhibitory 
edge from “DISC” to “MCLl”, respectively. These four nodes are indicated using solid black circles in 
Fig. [Ha). The corresponding plots on a logarithmic scale in the insets of (b) and (c) suggest a power-law 
scaling behavior between Tm and Ae [Eq. ©J. The fitted power-law scaling exponents axe /3 ^ 0.42, and 
0.38, respectively, for (b) and (c). 


perturbation exists that drives the system from the normal node to a caneerous node. 

Control implementation based on attractor network. Given a nonlinear dynamieal network 
in the real (physieal) space, the underlying phase space dimension may be quite high, rendering 
analysis of the dynamical behaviors difficult. The attractor network is a coarse grained represen¬ 
tation of the phase space, retaining information that is most relevant to the control task of driving 
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the network system to a desired final state. Once an attractor network has been constructed, actual 
control can be carried out through temporary changes in a set of experimentally adjustable param¬ 
eters, one at a time. This should be contrasted to one existing approach ll^ that requires accurate 
adjustments in the state variables, which may not always be realistic. 

We detail how actual control can be implemented based on the attractor network for the T-cell 
network. To be concrete, we assume that the control signal has the shape of a rectangular pulse in 
the plot of a parameter versus the time, as shown in Fig. [2a), where the control parameter is /r and 
the rectangular pulse has duration r and amplitude A/i = |/in — /^o|, with /io denoting the nominal 
parameter value and /i„ being the value during the time interval when control is on. For the T-cell 
network, we set /xq = 1.0. As /x is reduced the system approaches a bifurcation point. (In other 
examples a bifurcation can be reached by increasing a control parameter, as in low-dimensional 
GRNs detailed in Control Analysis.) Extensive numerical simulations in controlling the T-cell 
network from a cancerous state (Ci or C 2 ) to the normal state N shows that, to achieve control, 
there are wide ranges of choices for A/x and r. In fact, once /x„ is decreased through the bifurcation 
point He at which the initial attractor loses its stability, the goal of control can be realized. The 
critical value Hc for each parameter can be identified from a bifurcation analysis. Additionally, for 
Hn < /To there exists a required minimum control time r^, over which the system will move into 
the original basin of the target attractor before control is activated. Insofar as r > r^, the control 
signal can be released. Longer duration of control is not necessary since the system will evolve into 
the target attractor following its natural dynamical evolution associated with the nominal parameter 
/xq. The value of increases as /x„ is closer to Hc, where if /x„ = Hc, An is infinite due to the 
critical slowing down at the bifurcation point Hc- Figures I2b) and|2c) show, respectively, for the 
T-cell network, the relationship between and /x„ in controlling the strength of the activation 
edge from node “SIP” to node “PDGFR”, and that of the inhibitory edge from node “DISC” to 
node “MCLl” [cf.. Fig. [TJa), the nodes denoted as black circles and connected by bold coupling 
edges]. The critical value Hc (indicated by the dotted line) can be estimated accordingly. The insets 
of (b) and (c) show the corresponding plots of the relationships on a double logarithmic scale, with 
the horizontal axis to be Ae = /Xc — /x„, the exceeded value of /x„ over the critical point Hc- We 
observe the following power-law scaling behavior: 

An Q:|/^n dc\ y (2) 

where /3 is the scaling exponent. The region of control parameters at the upper-right region over 
the curve of rm(Ae), i.e., larger Ae value or longer duration r, corresponds to the case where 
control is successful in the sense that the system can definitely be driven to the desired final state. 

The power-law scaling relation for Tm demonstrated in Figs.jjb) and[2c) for the T-cell network 
is quite general, as it also holds for two-node and three-node GRNs (see Control Analysis). For 
the T-cell system, the critical values of parameters for all the possible controllable edges from 
Cl or C 2 to N, and the corresponding values of a and jd in Eq. © are provided in Table SI in 
Supplemental Information). The control magnitude and time for some parameters are identical, 
for the reason that the logic relationship from the corresponding edges to the same node can 
be described as “AND” (c.f.. Fig. [T]) so that in the continuous-time differential equation model, 
all these in-edges are equivalent. (The control results from the two-node and three-node GRNs 
between any pair of nearest-neighbor attractors are listed in Tables S2 and S3 in SI, respectively.) 
Due to the flexibility in choosing the control signal, our control scheme based on the attractor 
network is amenable to experimental implementation. 
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FIG. 3: Beneficial role of noise in controlling the T-cell network: probability that the normal state 
can be reached. Success rate to control the T-cell system from the cancerous state Ci to the normal state 
N using a combination of parameter perturbation and external noise (of amplitude) cr, where — Fc 

is the parameter deficiency. Warm colors indicate higher probability values of successful control. The 
perturbation duration is r = 200. The results are averaged over 1000 realizations. 


Beneficial role of noise in control. More than three deeades of intense researeh in nonlinear 
dynamieal systems has led to great knowledge about the role of noise, in terms of phenomena 
such as stochastic resonance llbOMbSll . coherence resonance feoMbOll . and noise-induced chaos OSll . 
etc. Common to all these phenomena is that a proper amount of noise can in fact be beneficial, 
for example, to optimize the signal-to-noise ratio, to enhance the signal regularity or temporal 
coherence, or to facilitate the transitions among the attractors. We find that, in our attractor- 
network based control framework, noise can also be beneficial. This can be understood intuitively 
by noting that our control mechanism is to make the system leave an undesired attractor and 
approach a desired one but noise in combination with parameter adjustments can facilitate the 
process of escaping from an attractor. To demonstrate this, we assume that the T-Cell network 
is subject to Gaussian noise, which can be modeled by adding independent normal distribution 
terms N(0, to the system equations, where a is the noise amplitude. We find that, with noise, 
the required magnitude of parameter change can be reduced. In fact, even when the controlled 


8 














parameter Hn has not yet reached the bifurcation point Hc, noise can lead to a non-zero probability 
for the system to escape the basin of the undesired attractor. 

Suppose the control parameter is set to the value /i„, which is insufficient to induce escape from 
the undesired attractor without noise. When noise is present, the system dynamics is stochastic. 
To characterize the control performance, we use a large number of independent realizations with 
the same initial condition. Specifically, we perform independent simulations starting from one 
cancerous state, e.g., Ci, but with insufficient control strength as characterized by the deficiency 
parameter Ad = /in — /io and calculate the probability P of control success through the number 
of trials that the system can be successfully driven to the normal state N. Figure [3] shows, 
on a double logarithmic scale, the values of P in the parameter plane of a and Ad, where the 
control parameter is the strength of the activation edge from node “SIP” to node “PDGFR” 
in the T-cell network. We see that, for fixed a, P decreases with Ad but, for any fixed value 
of Ad, the probability P increases with a, indicating the beneficial role of noise in facilitating 
control. In the parameter plane there exists a well-defined boundary, below which the control 
probability assumes large values but above which the probability is near zero. Testing alternative 
control parameters yields essentially the same results, due to the simplicity of the attractor network 
for the T-cell system and the multiple directed edges from each cancerous state to the normal state. 

Control Analysis 

In spite of the simplicity of its attractor network, the original T-cell network itself is still 
quite complicated from the point of view of nonlinear dynamical analysis. To have a better 
understanding of our control mechanism, we study GRNs of relatively low dimensions and carry 
out a detailed analysis of the associated attractor networks. 

Attractor network for a two-node GRN. We use a two-node GRN to understand the dynamical 
mechanism underlying the attractor network. As shown in Fig. |4l the network contains two auto¬ 
activation nodes (genes) and together they form a mutual inhibitory circuit. Such a topology was 
shown to be responsible for the regulation of blood stem cell differentiation ll^ . In addition, it 
is conceivable that such topologies can be constructed with tunable inputs using synthetic biology 
approaches 

In a typical experimental setting, four coupling parameters can be adjusted externally through 
the application of repressive or inductive drugs. To demonstrate attractor network and control im¬ 
plementation, we consider the parameter regime in which the system has four stable steady states 
(attractors) that correspond to four different cell states during cell development and differentiation. 
In particular, the dynamical network can be mathematically described as 


Xi 

= ai ■ 


+ bl- 

s’" 

— k ■ xi 

s^ + x^ 

S’" + X^ 

X2 

= 02 • 

^2 

+ (^2 ■ 

s’" 

-k- X 2 

s’" -f-x^ 

s’" -f xj" 


where the dynamical variables (xi, 0 : 2 ) characterize the protein abundances of the genes products, 
k denotes the degradation rate of each gene, and the tunable parameters oi, 02, hi, and 62 represent 
the strengths of auto or mutual regulations. In a GRN, the dynamical behaviors of inhibition and 
activation are captured by the Hill function: f{x) = x^/{x^ -f s"') for activation and f{x) = 
s”/(x” -f s”) for inhibition, where the parameter s characterizes half activation (or inhibition) 
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FIG. 4: A two-node GRN. Simplified model of the two-node GRN, where the arrowhead and bar-head 
edges represent activation and inhibition regulation, respectively. The sawtooth lines denote the strength of 
the tunable edge. 


concentration (for x = s, the output is 0.5), and n quantifies the correlation between the input and 
output eoneentrations, where a larger value of n eorresponds to a stronger inhibition or aetivation 
effeet. For any specific GRN, the values of both s and n can be determined experimentally. For 
simplicity, we assume the system to be symmetrie in that the inhibition and activation share the 
same Hill function (i.e., with the same parameters s and n). To have four attractors, we set the 
auto aetivation strengths, oi and 02 , to 1.0, and mutual inhibition strengths, bi and 62 , to 0.2. The 
value of the degradation rate fc is set to 1 . 1 , taking into aeeount the effects of protein degradation 
and eell volume expansion. 

Figure |5] shows a partieular process of eontrolling the system from an initial state, denoted as 
A, in whieh both xi and X 2 have low abundance, to a final state B where xi and X 2 have high 
and low abundanee, respectively. From the bifurcation diagram [Fig. I3a)] with respect to the 
control parameter oi, we see that, as oi is inereased from 1.0 to 1.4, in the lower braneh, the 
initial attractor A is destabilized through a saddle-node bifureation. The bifurcation-based eontrol 
proeess is shown in Figs. I3b-d), where panel (b) exhibits the phase spaee of the system prior to 
eontrol (oi = 1.0). When control is activated so that oi is set to Oi = 1.4, the original basin 
of attraetion of attraetor A merges into the basin of an intermediate attractor B', and the system 
originally in A starts to migrate towards the intermediate attractor B', as indicated by the arrowed 
trajeetory in panel (e). Control perturbation upon oi ean be withdrawn onee the state of the system 
enters the region belonging to the original basin of the target attraetor B, after whieh the system 
spontaneously evolves into B for oi = 1.0, as shown in Fig.|5];d). 

To obtain a global pieture of all possible control outcomes, we eonstruct the attractor network 
for the two-node GRN system, assuming that three eontrol parameters: Oi, 02 and b 2 , are available 


10 


(a) 



FIG. 5: Control of the two-node GRN. (a) Bifurcation diagram with respect to the control parameter ai, 
where the red and gray solid lines denote the stable and unstable steady states, respectively. In the two 
parallel cross-sections (with dashed line boundaries) for ai = a? and ai = o”, the yellow and brown dots 
represent the corresponding stable attractors, respectively. In (b-d), gray dashed lines represent the basin 
boundaries; black solid circles and green crosses denote attractors and unstable steady states, respectively, 
(b) For the initial parameter setting, oi = a^, the system is at a low concentration state A, and the target 
state is B. (c) By changing oi from to a”, attractor A is destabilized and its original basin is absorbed 
into that of the intermediate attractor B', so the system approaches B'. (d) When control perturbation upon 
is released, the landscape recovers to that in (b). Once the system has entered the basin of the target state 
B during the process in (c), it will evolve spontaneously towards B. Parameters in simulation are a? = 1.0, 
a” = 1.4, tQ = 0,ti = 23, and t 2 = 40. 

for control. The corresponding bifurcation diagrams are shown in Figs. [6ta-c), from which all 
saddle-node bifurcations can be identified for control design. When all the attractors are connected 
with directed and weighted edges through the control processes, i.e., when none of the attractor is 
isolated, we obtain an attractor network, as shown in Fig. |6td). Specifically, the edge weight can 
be assigned by taking into account the key characteristics of control such as the critical parameter 
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FIG. 6 : Construction of attractor network for two-node GRN. (a-c) Bifurcation diagrams of the system 
with respect to the coupling parameters oi, 02 and b 2 , respectively, where each bifurcation point can be ex¬ 
ploited to design control, (d) The corresponding attractor network, in which each directed edge corresponds 
to an elementary control that is designed to steer the system from the original attractor to the directed one. 
The solid and dashed edges, respectively, denote the positive and negative changes in the corresponding 
control parameters, (e) Sequential control signals required to drive the system from attractor A to attractor 
C through the path A ^ B —> C. In simulation, the original parameter values are a? = 1.0 and 62 = 0.2. 
We set a” = 1.4, followed by setting = 4.2, and the respective durations of the parameter perturbation 
are Ta = 23 and = 32. 


strength /ic and the power-law scaling behavior of the required minimum control time Tm (see 
Supplementary Table 2). From the attractor network, we can find all possible control paths for any 
given pair of original and desired states. 

From Fig. [b^d), we see that the two-node GRN system is fully controllable since any of the 
coexisting attractors is reachable by applying proper sequential controls upon the available param¬ 
eters. The concept of attractor network is appealing because it provides a clear control scenario 
to drive the system from any initial attractor to any desired attractor. In fact, the attractor net¬ 
work provides a blueprint that can be used to design a proper combination of parameter changes 
to induce the so-called synergistic or antagonistic effects ll70ll . For example, A is not directly 
connected with C, neither is B directly connected to D. However, the system can be steered from 
A to B through perturbation on oi, and then from B to C through parameter perturbation on 62 , as 
shown in Fig. [ 6 te). Another example to demonstrate the need of multiple parameter perturbation 
is to control the system from B to D. A viable control path is B —)■ C —> D, which can be realized 
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FIG. 7: Illustration of pseudo potential landscape. “Pseudo” potential E of the two-node GRN system 
(a) for ai = 1.0 (Ad « 0.3549), a = 0.05 and (b) for ai = 1.3 (Ad « 0.0549), cj = 0.05. Regions of 
warm and cold colors indicate the states with large and small pseudo energies, respectively. 


through perturbation on parameters ( 62 , Oi)- We also see that the two B —)■ A ^ D paths ean be 
realized through parameter changes in (oi, 02 ) and (oi, 62 ), respectively. 

When multiple control paths exist from an initial attractor to a final one, a practical issue is 
to identify an optimal path that is cost effective and robust. The concept of weighted-shortest 
path can be used to address this issue. Particularly, the weights of edges can be determined from 
experimental considerations such as the cost, limitation in drug dose, the control duration time, 
etc. 

Potential landscape and beneficial role of noise in nonlinear control. The role of noise in 
facilitating control of a nonlinear dynamical network can be understood using the concept of po¬ 
tential landscape ll35il7li[721] or Waddington landscape 117311 in systems bi ology, w hich essentially 
determines the biological paths for cell development and differentiation [I74l-l76n - the landscape 
metaphor. The potential landscape has been used to manipulate time scales to control stochastic 
and induced switching in biophysical networks [f^ . Intuitively, the power of the concept of the 
landscape can be understood by resorting to the elementary physical picture of a ball moving in 
a valley under gravity. The valley thus corresponds to one stable attractor. To the right of the 
valley there is a hill, or a potential barrier in the language of classical mechanics. The downhill 
side to the right of the barrier corresponds to a different attractor. Suppose the confinement of 
ball’s motion within the valley is undesired and one wishes to push the ball over the barrier to the 
right attractor (desired). If the valley is deep (or the height of the barrier is large), there will be 
little probability for the ball to move across the top of the barrier towards the desired attractor. In 
this case, a small amount of noise is unable to enhance the crossover probability. However, if the 
barrier height is small, a small amount of noise can push the ball over to desired attractor on the 
right side of the barrier. Thus, the beneficial role of noise is more pronounced for small height 
of the potential barrier, a behavior that we observe when controlling the T-cell network (Fig. (3]). 
In mechanics, the system can be formulated using a potential function so that, mathematically. 
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the motion of the ball ean be described by the Langevin equation, which has been a paradigmatic 
model to understand nonlinear phenomena such as stochastic resonance [6^^]. In the past few 
years, a quantitative approach has been developed to mapping out the potential landscape for gene 


circuits or gene regulatory networks [l35Ll77l-l79ll. In nonlinear dynamical systems, a similar con¬ 


cept exists - quasipotential |8^82], which plays an important role in understanding phenomena 
such as noise-induced chaos. 

For an attractor network, in the presence of noise each node corresponds to a potential valley 
of certain depth characterizing the stability of the attractor. For a fixed depth, noise of larger 
amplitude a leads a larger escaping probability or shorter escaping time. When the amplitude of 
the control signal is not sufficient to drive the system across the local potential barrier, noise can 
facilitate control by pushing the system out of the undesired valley (attractor). 

The potential landscape for a GRN under Gaussian noise can be constructed from the dynam¬ 
ical equations of the system using the concept of “pseudo” energy iIt^ (see Methods). For the 
two-node GRN system [Eq. ([3])] subject to stochastic disturbance iV(0, cr^), we can compute the 
potential landscape for any combination of a system parameter (say Oi) and the noise amplitude 
a. Figure 0 shows two examples of the landscapes for oi = 1.0 and Oi = 1.3, where the noise 
amplitude is a = 0.05. We see that, for oi = 1.0, there are four valleys (attractors). For oi = 1.3, 
the pseudo energy for A (the original valley at the lower-left comer) becomes higher, and the 
path for the transition from A to B becomes more pronounced. Further increasing oi towards the 
critical value (about 1.35) raises the energy of A to the level of the potential barrier, effectively 
eliminating the corresponding valley and the attractor itself. 

Attractor network for a three-node GRN. We also study a three-node GRN system, as shown 
in Fig. [8]; a). Similar to the two-node GRN system, there exist both auto and mutual regulations 
among the nodes. All the interactions are assumed to be characterized by the same parameters of 
s and n in the Hill function. The nonlinear dynamical equations of the system are 11471 18311 
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where the state variables (xi, X 2 and xf) represent the abundances of the three genes products, 
the auto-activation parameters Oi, b 2 , C3 and the mutual-inhibition parameters 02 , 03, bi, 63, Ci, 
C 2 are all experimentally accessible. To be concrete, initially all the auto activation and mutual 
inhibition parameters are set to be 1.0 and 0.1, respectively, and k is the degradation rate that can 
be conveniently set to unity. The parameters in the Hill function are n = 4 and s = 0.5. There are 
altogether eight attractors in this system, as shown in Fig.[8tb), which are distributed symmetrically 
in the three-dimensional state space. For example, attractor H has relatively high values for all 
three dynamical variables, and attractor B exhibits the opposite case with low abundances. For 
attractors A, C and F, one of the three state variables is high and the other two are low. For 
attractors D, E and G, one of the three state variables is low but the other two are high. 

From numerical simulations, we find that the features of control are essentially the same as 
those for the two-node GRN system, in terms of characteristics such as the existence of critical 
control strength and the power-law behavior of the minimum control time (see Table S3 in SI). 
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FIG. 8: Attractor network and control of a three-node GRN system, (a) Schematic illustration of a 
three-node GRN system. The arrowhead and bar-head edges represent activation and inhibitory regulations, 
respectively. The sawtooth lines specify that corresponding edge strength is experimentally adjustable, (b) 
Coexisting attractors (A to H) in the phase space, (c) The underlying attractor network, where each node 
represents an attractor and each weighted directed link indicates that its strength can be experimentally 
tuned to steer the system from the starting attractor to the pointed attractor. 


We construct the attractor network in Fig. [8tc) through combinations of all eight attractors (as 
nodes) and directed elementary controls (as weighted directed edges). Information in Table S3 
can also be used to estimate the respective weights of the edges. From the attractor network, for 
any given pair of initial and final states, we can identify all the viable control paths. Furthermore, 
the weighted-shortest path can be calculated once the edge weights are determined. 

We note that, typically, the attractor network based on elementary control is not an all-to-all 
directed network so that certain control paths are absent, e.g., from attractor H to B. The 
biological meaning is that, while a stem state can be differentiated into various types of cells 
through bifurcation, the opposite paths back to the stem state are much more difficult to find. 

Discussions 


The field of controlling chaos in nonlinear dynamical systems has been active for more 
than two decades since the seminal work of Ott, Grebogi, and Yorke [l84n . The basic idea was 
that chaos, while signifying random or irregular behavior, possesses an intrinsically sensitive 
dependence on initial conditions, which can be exploited for controlling the system using only 
small perturbations. This feature, in combination with the fact that a chaotic system possesses 
an infinite set of unstable periodic orbits, each leading to different system performance, implies 
that a chaotic system can be stabilized about some desired state with optimal performance using 
small control perturbations. Controlling chaos has since been studied extensively and examples 
of successful experimental implementation abound in physical, chemical, biological, and engi¬ 
neering systems llSSn . The vast literature on controlling chaos, however, has been mostly limited 
to low-dimensional systems, systems that possess one or a very few unstable directions (i.e., 
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one or a very few positive Lyapunov exponents). Complex networks with nonlinear dynamies 
are generally high dimensional, rendering inapplieable existing methodologies of ehaos eontrol. 
While mathematieal frameworks of eontrollability of eomplex networks [0, @] were developed 
and extensively studied reeently, the defieieney of sueh rigorous mathematieal frameworks is that 
the nodal dynamieal proeesses must be assumed to be linear. While eontrollability for nonlinear 
eontrol ean be formulated based on Lie braekets 11861] . it may be diffieult to implement the abstraet 
framework for eomplex networks. 

Controlling nonlinear dynamies on eomplex networks remains to be an outstanding and ehal- 
lenging problem, espeeially in terms of the two key issues: eontrollability and aetual eontrol. To 
assess the eontrollability of nonlinear dynamieal networks, drastieally different approaehes than 
the linear eontrollability framework are needed. While there were previous works on speeifie eon- 
trol methods sueh as pinning eontrol [2^22] and brute-foree eontrol that rely on altering the state 
variables of the underlying system (whieh in realistie situations ean be diffieult to implement), we 
eontinue to laek a general framework for aetual eontrol of eomplex networks with nonlinear dy¬ 
namies through realistie, physieal means. The main diffieulty lies in the extremely diverse nonlin¬ 
ear dynamieal behaviors that a network ean generate, making it praetieally impossible to develop a 
general mathematieal framework for eontrol. In partieular, the traditional eontrol theoretieal tools 
for linear dynamieal systems aim to eontrol the detailed states of all of the variables, whieh is 
in faet an overkill for most systems. For nonlinear dynamieal networks, a physieally meaningful 
approaeh may not require detailed eontrol of all state variables. With this relaxation of the eontrol 
requirement, it may be possible to develop a framework of eontrollability and devise aetual eontrol 
strategies for nonlinear dynamieal networks based on physieal/experimental eonsiderations. 

A eommon feature of nonlinear dynamieal systems is the emergenee of a large number of 
distinet, eoexisting attraetors 11281.18711 . Often the performanee and funetions of the system are 
determined by the partieular attraetor that the system has settled into, assoeiated with which the 
detailed states of the dynamical variables are not relevant. The key is thus to develop control 
principles whereby we nudge a complex, nonlinear system from attractor to attractor through small 
perturbations to a set of physically or experimentally feasible parameters. The main message of 
this paper is that a controllability framework can be developed for nonlinear dynamical networks 
based on the control of attractors. 

Generally, the reason for control is that the current system is likely to evolve into an undesired 
state (attractor) or the system is already in such a state, and one wishes to apply perturbations 
to bring the system out of the undesired state and steer it into a desired state. The first step is 
then to identify a final state or attractor of the system that leads to the desired performance. The 
next step is to choose a set of experimentally adjustable parameters and determine whether small 
perturbations to these parameter can bring the system to the desired attractor. That is, under physi¬ 
cally realizable perturbations there should be a control path between the undesired and the desired 
attractors. The path can be directly from the former to the latter, or there can be intermediate at¬ 
tractors on the path. For example, due to the physical constraint on the control parameters and the 
ranges in which they can be meaningfully varied, one can drive the system into some intermediate 
attractors by perturbing one set of parameters, and then from these attractors to the final attractor 
by using a different set of parameter control. For a complex, nonlinear dynamical network, the 
number of coexisting attractors can be large. Given a set of system performance indicators, one 
can classify all the available attractors into three categories: the undesired, desired, and the inter¬ 
mediate attractors. We say a nonlinear network is controllable if there is a control path from any 
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undesired attractor to the desired attractor under finite parameter perturbations. Regarding each 
attractor as a node and the control paths as directed links or edges, we generate an attractor net¬ 
work whose properties determine the controllability of the original networked dynamical system. 
For example, the average path length from an undesired to a desired attractor and the “control 
energy” (or the amount of necessary parameter perturbations) can serve as quantitative measures 
to characterize the controllability of the original network. We demonstrate our idea of control and 
construction of attractor networks using realistic networks from systems and synthetic biology. 
We also find that noise can facilitate control of nonlinear dynamical networks, and we provide a 
physical theory to understand this counterintuitive phenomenon. 

While we emphasize the need to focus on physically meaningful and experimentally accessible 
parameter perturbations, there can still be a large number of attractor networks depending on the 
choice of the parameters, making it difficult to formulate a rigorous mathematical framework. 
We believe that these issues can and will be satisfactorily addressed in the near future, ultimately 
realizing the grand goal of controlling nonlinear dynamical networks. 

Methods 


Pseudo potential landscape. For a dissipative, nonlinear dynamical system subject to noise, 
we can construct a pseudo potential landscape based on the state probability distribution. Assume 
that, asymptotically, the system approaches a stationary distribution. For a canonical dynamical 
system, the potential can be defined as E{x) = — log F’(x), where F’(x) is the probability density 
function. For a conservative dynamical system, the direction of system evolution is nothing but 
the direction of the gradient of the potential function. However, this does not hold for dissipative 
dynamical systems. The potential function thus does not have the same physical meaning as that 
for a conservative system, henceforth the term pseudo potential. This approach can be adopted to 
GRNs. 

To obtain the stationary distribution, we use the modified weighted-ensemble algorithm pro¬ 
posed by Kromer et al. 11721] . which offers faster convergence than, for example, the traditional 
random walk method. To be illustrative, we take the two-node GRN system [Eq. ([3])] as an ex¬ 
ample to demonstrate how the pseudo potential landscape can be numerically obtained. The state 
space of the two-dimensional dynamical system is partitioned into an M x M lattice with reflective 
boundaries conditions. Initially the probability Pm,nil) of all gird points are set to be uniform. The 
simulation time is divided into T steps, where each step has the duration r. At the beginning of 
each step t, there are N walkers randomly distributed at the grid point (m, n), which carry equal 
weight Pm,nil)/N and perform random walk under the system dynamics and noise. The locations 
of these walkers in the grid are recorded at the end of each time step, and the probability at next 
time step, Pm,nil + 1), is the summation of the probabilities carried by all the walkers at time t. At 
time it + 1), N new walkers carrying the updated probability at each grid point perform random 
walk again on the grid. This procedure repeats until the probability distribution becomes station¬ 
ary, say Pm,n, which gives the pseudo potential landscape as Eim, n) = — log Pm,n- Numerically, 
the time evolution of all walkers can be simulated using the second-order Heun method for inte¬ 
grating stochastic differential equations. For Fig.|71 the state space is divided into a 500 x 500 grid. 
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At each grid point there are = 20 walkers, each evolving T = 2000 time steps with r = 10 
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